;; Non-robust row reduced echelon form functions
;; Copyright (C) 2025  Rocks Mazama
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <https://www.gnu.org/licenses/>.

(ns mazama.linalg.rref
  (:require
   [mazama.linalg.r-e-f :as r-e-f]))

;; See mazama.linalg.matrix for a description of the map used to describe
;; matrices.
;;
;; Augmented matrices can also be used. Augmented matrices of the form [A|b]
;; have a key :acols the value of which is the number of columns in the matrix
;; A. :ncols represents the number of columns in the entire augmented matrix.

(defn- rref-step
  [rows pivot-column-index ero-scale]
  (let [pivot-row (first rows)
        pivot-value (nth pivot-row pivot-column-index)
        f (fn [i row]
            (if (zero? i)
              (map #(/ % pivot-value) row)
              (let [v (nth row pivot-column-index)
                    sf (/ v pivot-value)]
                (map #(- %1 (* sf %2)) row pivot-row))))]
    {:pivot-rows (map-indexed f rows) :ero-scale (* ero-scale pivot-value)}))

(defn rref
  "Compute the reduced row echelon form of the matrix `m`.
  `m` can be an augmented matrix (with key :acols)."
  [m]
  (let [{:keys [pivot-column-indexes pivot-rows zero-rows ero-scale]} (r-e-f/r-e-f m ())]
    (loop [pivot-rows (seq pivot-rows)
           results ()
           pci pivot-column-indexes
           ero-scale ero-scale]
      (if pivot-rows
        (let [{:keys [pivot-rows ero-scale]} (rref-step pivot-rows (first pci) ero-scale)]
          (recur (next pivot-rows)
                 (conj results (first pivot-rows))
                 (next pci)
                 ero-scale))
        {:pivot-column-indexes (reverse pivot-column-indexes)
         :pivot-rows results
         :zero-rows zero-rows
         :ero-scale ero-scale}))))
